!---------------------------
! EXPLANATION OF THE MODULE
!---------------------------
! THIS MODULE CONTAINS INTERFACES FOR ALL EXTERNAL SUBROUTINES AND FUNCTIONS
! THIS MODULE IS MEANT FOR FUTURE EXPANSION AS AN EMPTY SLOT. ANY NEW SUBROUTINES
! OR FUNCTIONS OF GENERAL PURPOSE SHOULD BE WRITTEN HERE

MODULE ALL_INTERFACES_2D

    IMPLICIT NONE

    INTEGER,PARAMETER::DBL=SELECTED_REAL_KIND(15,307)      
    INTEGER,PARAMETER::SIZE_PROB=2                                    ! SIZE OF PROBLEM: 2D OR 3D
    INTEGER,PARAMETER::NSURFACE=35                                    ! NUMBER OF SURFACES
    INTEGER,PARAMETER::NDIVISION=100                                ! CONJUGATE POINT SOLUTION
    REAL(KIND=DBL),PARAMETER::TOLERANCE=EPSILON(1.0D0)*100.0D0

    INTEGER,PARAMETER::MAXSIZE=2

    TYPE,PUBLIC::TENSOR_2R
        !PRIVATE
        INTEGER::TSIZE
        REAL(KIND=DBL),DIMENSION(MAXSIZE,MAXSIZE)::COMPONENT
    END TYPE TENSOR_2R

    TYPE,PUBLIC::TENSOR_4R
        INTEGER::TSIZE
        REAL(KIND=DBL),DIMENSION(MAXSIZE,MAXSIZE,MAXSIZE,MAXSIZE)::COMPONENT
    END TYPE TENSOR_4R

    ! STRESS AND STRAIN   
    TYPE(TENSOR_2R)::STRESS_TENSOR_CURRENT, STRESS_TENSOR_NEXT, STRAIN_TENSOR_CURRENT, STRAIN_TENSOR_NEXT, &
    & PLASTIC_STRAIN_INCREMENT
    
    ! TANGENTS
    TYPE(TENSOR_4R)::CONSTITUTIVE_TENSOR
    REAL(KIND=DBL)::EL_TANGENT(3,3),ELPL_TANGENT_CURRENT(3,3),ELPL_TANGENT_NEXT(3,3)

    ! MODEL PARAMETERS
    REAL(KIND=DBL)::VOID_RATIO, PATM, PMIN, PREF
    REAL(KIND=DBL)::SHEAR_CONST_A, SHEAR_CONST_N, POISSON
    REAL(KIND=DBL)::ETAMAX_CONST_A1, ETAMAX_CONST_B1
    REAL(KIND=DBL)::GNMAX_CONST_A2, GNMAX_CONST_B2
    REAL(KIND=DBL)::GNMIN_CONST_A3, GNMIN_CONST_B3,DEGRADATION
    REAL(KIND=DBL)::DILATANCY_COEFF_MU_MONO, DILATANCY_COEFF_MU_CYCLIC, DILATANCY_STRAIN_SC, CRITICAL_STRESS_RATIO_M
    REAL(KIND=DBL)::INITIAL_STRAIN

    ! COORDINATES OF THE STEADY STATE LINE IN VOID RATIO - HYDROSTATIC PRESSURE SPACE
    REAL(KIND=DBL)::STEADY_STATE_DATA(10,2)

    ! COORDINATES OF THE HYDROSTATIC COMPRESSION LINE IN VOID RATIO-HYDROSTATIC PRESSURE SPACE
    REAL(KIND=DBL)::HYDROSTATIC_DATA(10,2)

    ! PARAMETERS OF THE HARDENING SURFACES
       TYPE(TENSOR_2R),DIMENSION(NSURFACE)::CENTRE,REVERSAL
    REAL(KIND=DBL)::GAMMA(NSURFACE),MAX_STRAIN,GEFFECTIVE
    INTEGER::SURFACE,KEY

    ! PLASTICITY PARAMETERS
    REAL(KIND=DBL)::HYDROSTATIC_PRESSURE, ELASTIC_SHEAR_MODULUS, GNMAX, GNMIN, ETAMAX, LAMDA, RADIUS, &
    & PLASTIC_SHEAR_MODULUS, PLASTIC_MODULUS_REVERSAL,DILATANCY,DILATANCY_AT_REVERSAL

    ! STATE PARAMETERS
    REAL(KIND=DBL),PRIVATE::STATE_INDEX, STEADY_STATE_VOID_RATIO, HYDROSTATIC_VOID_RATIO, DELTA_VOLUMETRIC_STRAIN

    INTERFACE RECEIVEDATA
        MODULE PROCEDURE RECEIVEDATA_REAL,RECEIVEDATA_INT,RECEIVEDATA_REAL_A
    END INTERFACE

    INTERFACE SENDDATA
        MODULE PROCEDURE SENDDATA_REAL,SENDDATA_INT,SENDDATA_REAL_A
    END INTERFACE

    INTERFACE SOLVECONJUGATE
        MODULE PROCEDURE SOLVECONJUGATE_TYPE1
    END INTERFACE
 
    ! TENSOR OPERATORS
    INTERFACE OPERATOR (+)
        MODULE PROCEDURE TSUM
    END INTERFACE
    
    INTERFACE OPERATOR (-)
        MODULE PROCEDURE SUBTRACT                            
    END INTERFACE
    
    INTERFACE OPERATOR (.DDOT.)
        MODULE PROCEDURE CONTRACT                            
    END INTERFACE
    
    INTERFACE OPERATOR (*)
        MODULE PROCEDURE MULTIPLY_DOUBLE_REAL, MULTIPLY_SINGLE_REAL                        
    END INTERFACE

    INTERFACE TENSORMULTIPLY
        MODULE PROCEDURE TENSORMULTIPLY222, TENSORMULTIPLY422
    END INTERFACE

        ! EXPLANATION OF THE VARIABLES:
        
        ! STRESS_CURRENT        :    ARRAY OF 6/3 DIMENSIONS TO STORE CURRENT VALUES OF CURRENT STRESS
        ! STRAIN_CURRENT         :    ARRAY OF 6/3 DIMENSIONS TO STORE CURRENT VALUES OF CURRENT STRAIN
        ! STRAIN_NEXT             :    ARRAY OF 6/3 DIMENSIONS TO STORE CURRENT VALUES OF NEXT STRAIN
        
        ! MODEL_PARAMETERS AS PER: CUBRINOVSKI AND ISHIHARA (1998)
        
        ! MODEL_PARAMETER(1)    :    VOID_RATIO
        ! MODEL_PARAMETER(2)    :    SHEAR_CONST_A
        ! MODEL_PARAMETER(3)    :    SHEAR_CONST_N
        ! MODEL_PARAMETER(4)    :    POISSON
        ! MODEL_PARAMETER(5)    :    ETAMAX_CONST_A1
        ! MODEL_PARAMETER(6)    :    ETAMAX_CONST_B1
        ! MODEL_PARAMETER(7)    :    GNMAX_CONST_A2
        ! MODEL_PARAMETER(8)    :    GNMAX_CONST_B2
        ! MODEL_PARAMETER(9)    :    GNMIN_CONST_A3
        ! MODEL_PARAMETER(10)   :    GNMIN_CONST_B3
        ! MODEL_PARAMETER(11)   :    DEGRADATION
        ! MODEL_PARAMETER(12)   :    DILATANCY_COEFF_MU_MONO
        ! MODEL_PARAMETER(13)   :    DILATANCY_COEFF_MU_CYCLIC     
        ! MODEL_PARAMETER(14)   :    DILATANCY_STRAIN_SC
        ! MODEL_PARAMETER(15)   :    CRITICAL_STRESS_RATIO_M
        ! MODEL_PARAMETER(16)   :    PATM

        ! SSL_VOID_RATIO(10)    :    STEADY STATE LINE DATA - VOID RATIO VERSUS EFFECTIVE PRESSURE
        ! SSL_PRESSURE(10)
        
        ! HSL_VOID_RATIO(10)    :    HYDROSTATIC LINE DATA - VOID RATIO VERSUS EFFECTIVE PRESSURE 
        ! HSL_PRESSURE(10)            
        
        ! HARDENING_PARAMETER_INT(1)                             : SURFACE
        ! HARDENING_PARAMETER_INT(2)                             : KEY
        
        ! HARDENING_PARAMETER_REAL(1:3*NSURFACE)                 : CENTRE
        ! HARDENING_PARAMETER_REAL(3*NSURFACE+1:6*NSURFACE)      : REVERSAL_ORDINATE
        ! HARDENING_PARAMETER_REAL(6*NSURFACE+1:7*NSURFACE)      : PLASTIC_STRAIN
        ! HARDENING_PARAMETER_REAL(7*NSURFACE+1)                 : MAX_STRAIN
        ! HARDENING_PARAMETER_REAL(7*NSURFACE+2)                 : GEFFECTIVE
        ! HARDENING_PARAMETER_REAL(7*NSURFACE+3)                 : INITIAL_STRAIN
     
    CONTAINS

    SUBROUTINE MODEL_2D (STRESS_CURRENT, STRAIN_CURRENT, STRAIN_NEXT, MODEL_PARAMETER, SSL_VOID_RATIO, &
    & SSL_PRESSURE, HSL_VOID_RATIO, HSL_PRESSURE, HARDENING_PARAMETER_REAL, HARDENING_PARAMETER_INT,TANGENT)                
    
    REAL(KIND=DBL),DIMENSION(3),INTENT(IN)::STRAIN_NEXT
    REAL(KIND=DBL),DIMENSION(10),INTENT(IN)::SSL_VOID_RATIO, SSL_PRESSURE
    REAL(KIND=DBL),DIMENSION(10),INTENT(IN)::HSL_VOID_RATIO, HSL_PRESSURE
    REAL(KIND=DBL),DIMENSION(3),INTENT(INOUT)::STRAIN_CURRENT
    REAL(KIND=DBL),DIMENSION(3),INTENT(INOUT)::STRESS_CURRENT
    REAL(KIND=DBL),DIMENSION(16),INTENT(INOUT)::MODEL_PARAMETER
    REAL(KIND=DBL),DIMENSION(7*NSURFACE+3),INTENT(INOUT)::HARDENING_PARAMETER_REAL
    INTEGER,DIMENSION(2),INTENT(INOUT)::HARDENING_PARAMETER_INT
    REAL(KIND=DBL),DIMENSION(3,3),INTENT(INOUT)::TANGENT
    REAL(KIND=DBL)::PRESSURE,CURRENT_STRAIN,CURRENT_RADIUS,CURRENT_MU0,GN,K1_HP,K2_MU,TEMP(2,2)
    TYPE(TENSOR_2R)::SIJ, DEPSILON, GIJ, FIJ, SIJ_NEXT, NIJ, XIJ_CP
    LOGICAL::ERROR ! IN DLL VERSION
      
    !********************************************************************************************
    ! ASSIGN TO GLOBAL VARIABLES
    !********************************************************************************************
    
    CALL RECEIVEDATA(DATAIN=STRESS_CURRENT,DSIZE=3,CHOICE=1)     
    CALL RECEIVEDATA(DATAIN=STRAIN_CURRENT,DSIZE=3,CHOICE=2)
    CALL RECEIVEDATA(DATAIN=STRAIN_NEXT,DSIZE=3,CHOICE=3)
    CALL RECEIVEDATA(DATAIN=MODEL_PARAMETER,DSIZE=16,CHOICE=4)
    CALL RECEIVEDATA(DATAIN=SSL_VOID_RATIO,DSIZE=10,CHOICE=5)
    CALL RECEIVEDATA(DATAIN=SSL_PRESSURE,DSIZE=10,CHOICE=6)
    CALL RECEIVEDATA(DATAIN=HSL_VOID_RATIO,DSIZE=10,CHOICE=7)
    CALL RECEIVEDATA(DATAIN=HSL_PRESSURE,DSIZE=10,CHOICE=8)
    CALL RECEIVEDATA(DATAIN=HARDENING_PARAMETER_REAL,DSIZE=7*NSURFACE+3,CHOICE=9)
    CALL RECEIVEDATA(DATAIN=HARDENING_PARAMETER_INT,DSIZE=2)
    CALL RECEIVEDATA(DATAIN=TANGENT,DSIZE=3)

    !*********************************************************************************************
    ! NEGATIVE EFFECTIVE PRESSURE - ABORT
    !*********************************************************************************************

    PRESSURE=FIRSTINVARIANT(STRESS_TENSOR_CURRENT)
    IF(PRESSURE.LE.TOLERANCE) THEN
          TEMP = 0D0
          STRESS_TENSOR_NEXT = CREATETENSOR(TEMP,SIZE_PROB)
          ELPL_TANGENT_NEXT = 0D0
          PRINT *,'STRESS DILATANCY MODEL::ALL INTERFACES::NEGATIVE HYDROSTATIC PRESSURE'
          CALL UPDATE_PARAMETERS()
        RETURN
    END IF

    !*********************************************************************************************
    ! NO CHANGE IN STRAIN - RETURN
    !*********************************************************************************************
    
    DEPSILON=STRAIN_TENSOR_NEXT-STRAIN_TENSOR_CURRENT
    IF (NORM(DEPSILON).LE.TOLERANCE) THEN
        WRITE(*,*)'STRESS DILATANCY MODEL::ALL INTERFACES, NO CHANGE IN STRAIN'
        RETURN
    END IF
   
    !***********************************************************************************************
    ! CALCULATE GE, GNMAX, GNMIN, ETAMAX AND RANK 4 CONSTITUTIVE TENSOR
    !***********************************************************************************************
    
    CALL GET_REF_PRESSURES()
    CALL CALCULATE_ELASTIC_SHEAR_MODULUS()
    CALL CALCULATE_GNMAX()
    CALL CALCULATE_GNMIN()
    CALL CALCULATE_ETAMAX()
    CALL CALCULATE_RADIUS()
    CALL SETINITIALTANGENT()
    CONSTITUTIVE_TENSOR=CONSTITUTIVETENSOR(ELASTIC_SHEAR_MODULUS,POISSON,SIZE_PROB)
    
    !***********************************************************************************************
    ! CHECK IF THE STRAIN INCREMENT IS HYDROSTATIC IN NATURE
    !***********************************************************************************************
    
    !DEPSILON = STRAIN_TENSOR_NEXT-STRAIN_TENSOR_CURRENT
    !IF(NORM(DEVIATORIC(DEPSILON)).LE.TOLERANCE) THEN
    !    WRITE(*,*)'STRESS DILATANCY MODEL::ALL INTERFACES, HYDROSTATIC STRAIN INCREMENT'
    !    STRESS_TENSOR_NEXT=STRAIN_TENSOR_CURRENT+TENSORMULTIPLY(CONSTITUTIVE_TENSOR, DEPSILON)
    !    ELPL_TANGENT_NEXT = EL_TANGENT
    !    CALL UPDATE_PARAMETERS()
    !    RETURN
    !END IF
        
    !***********************************************************************************************
    ! PROJECT EXISTING STRESS TO A REFERENCE PLANE
    !***********************************************************************************************
    
    SIJ=DEVIATORIC(STRESS_TENSOR_CURRENT)
    SIJ=PREF/HYDROSTATIC_PRESSURE*SIJ

    !**********************************************************************************************
    ! INITIALISE PLASTIC MODULUS, DILATANCY & YIELD SURFACE AT THE START OF LOADING WHEN KEY = 0
    !**********************************************************************************************

    IF (KEY.EQ.0) THEN

        !WRITE(*,*),'KEY = ',KEY
        SURFACE                  = 1
        PLASTIC_SHEAR_MODULUS    = GNMAX*HYDROSTATIC_PRESSURE
        PLASTIC_MODULUS_REVERSAL = PLASTIC_SHEAR_MODULUS
        DILATANCY                = DILATANCY_COEFF_MU_MONO        
        REVERSAL(SURFACE)        = SIJ
        
        ! GET ANISOTROPIC STRAIN FROM THE BACKBONE CURVE
        INITIAL_STRAIN = CALCULATE_EQV_STRAIN(SECINVDEV(SIJ))
        !WRITE(*,*)'Q',SECINVDEV(SIJ)
        !WRITE(*,*)'P',HYDROSTATIC_PRESSURE
        !WRITE(*,*)'INITIAL_STRAIN',INITIAL_STRAIN
        !INITIAL_STRAIN = CURRENT_STRAIN
        
        ! NOW ASSIGN THIS TO ALL LOADING SURFACES AS THE INITIAL STRAIN
        GAMMA(1:NSURFACE)=INITIAL_STRAIN        
                     
        ! LOCATE CENTRE OF LOADING SURFACE 1
        NIJ = NORMALISE(SIJ)
        XIJ_CP=SOLVECONJUGATE(SIJ,NIJ)
        CENTRE(SURFACE)=SIJ-(NORM(SIJ)/NORM(XIJ_CP))*(XIJ_CP)
        !WRITE(*,*)'SIJ',SIJ
        !WRITE(*,*)'XIJ_CP',XIJ_CP
        !WRITE(*,*)'NORMSIJ',NORM(SIJ)
        !WRITE(*,*)'NORMXIJ_CP',NORM(XIJ_CP)
        !WRITE(*,*)'CENTRE',CENTRE(SURFACE)
      
        ! CALCULATE TANGENT
        ELPL_TANGENT_CURRENT = EL_TANGENT

        ! SET CURRENT DILATANCY COEFF AS MU_MONO
        CURRENT_MU0 = DILATANCY_COEFF_MU_MONO

        DILATANCY_AT_REVERSAL=CURRENT_MU0+(CRITICAL_STRESS_RATIO_M-CURRENT_MU0)* &
        & REAL(ATAN(INITIAL_STRAIN/DILATANCY_STRAIN_SC),KIND=DBL)/REAL(ASIN(1.0D0),KIND=DBL)
   
    ELSE

       !WRITE(*,*),'KEY = ',KEY
       !WRITE(*,*),'SURFACE = ',SURFACE

       CURRENT_RADIUS = NORM(REVERSAL(SURFACE)-CENTRE(SURFACE))  
       CURRENT_STRAIN = GAMMA(SURFACE)
    
       IF(SURFACE.EQ.1) THEN
           K1_HP = 1.0D0
           K2_MU = 1.0D0
           CURRENT_MU0 = DILATANCY_COEFF_MU_MONO
       ELSE
           K1_HP = 2.0D0
           IF (CURRENT_STRAIN.LE.0.001D0) THEN
           !IF (CURRENT_STRAIN.LE.0.05D0) THEN
               K2_MU = 2.0D0
           ELSEIF (CURRENT_STRAIN.GE.0.003D0) THEN
           !ELSEIF (CURRENT_STRAIN.GE.0.01D0) THEN
               K2_MU = 1.0D0
           ELSE 
               K2_MU = 2.0D0-(CURRENT_STRAIN-0.001D0)/0.002D0
               !K2_MU = 2.0D0-(CURRENT_STRAIN-0.05D0)/(0.1D0-0.05D0)
           END IF
           CURRENT_MU0 = DILATANCY_COEFF_MU_CYCLIC
       END IF

   !ENDIF
    
       PRESSURE = FIRSTINVARIANT(STRESS_TENSOR_CURRENT)/REAL(SIZE_PROB,KIND=DBL)
                
       GN=(GNMAX-GNMIN)*DEXP(-DEGRADATION*CURRENT_STRAIN/0.01D0/K1_HP)+GNMIN
   
       PLASTIC_SHEAR_MODULUS = (GN-(GN-GNMIN)*(DEGRADATION*CURRENT_STRAIN/0.01D0/K1_HP))* &
       & (1.0D0-CURRENT_RADIUS/RADIUS)**2*PRESSURE
   
       GN=(GNMAX-GNMIN)*DEXP(-DEGRADATION*INITIAL_STRAIN/0.01D0/K1_HP)+GNMIN 
   
       PLASTIC_MODULUS_REVERSAL = (GN-(GN-GNMIN)*(DEGRADATION*INITIAL_STRAIN/0.01D0/K1_HP))*PRESSURE
             
       DILATANCY=CURRENT_MU0+(CRITICAL_STRESS_RATIO_M-CURRENT_MU0)* &
       & REAL(ATAN(CURRENT_STRAIN/DILATANCY_STRAIN_SC/K2_MU),KIND=DBL)/REAL(ASIN(1.0D0),KIND=DBL)
            
       DILATANCY_AT_REVERSAL = DILATANCY_COEFF_MU_CYCLIC+(CRITICAL_STRESS_RATIO_M-DILATANCY_COEFF_MU_CYCLIC)* &
       & REAL(ATAN(INITIAL_STRAIN/DILATANCY_STRAIN_SC/K2_MU),KIND=DBL)/REAL(ASIN(1.0D0),KIND=DBL)  

       GN=(GNMAX-GNMIN)*DEXP(-DEGRADATION*CURRENT_STRAIN/0.01D0/K1_HP)+GNMIN

    ENDIF

    !WRITE(*,*)'PLASTIC MODULUS',PLASTIC_SHEAR_MODULUS
    !WRITE(*,*)'PLASTIC MODULUS AT REVERSAL',PLASTIC_MODULUS_REVERSAL
    !WRITE(*,*)'GN',GN
    !WRITE(*,*)'GNMAX',GNMAX
    !WRITE(*,*)'GNMIN',GNMIN
    !WRITE(*,*)'DILATANCY',DILATANCY
    !WRITE(*,*)'DILATANCY AT REVERSAL',DILATANCY_AT_REVERSAL
    !WRITE(*,*)'MU NOT',CURRENT_MU0
    !WRITE(*,*)'K2_MU',K2_MU
    !WRITE(*,*)'CURRENT_RADIUS',CURRENT_RADIUS
    !WRITE(*,*)'CURRENT_STRAIN',CURRENT_STRAIN
  
    CALL OPTIMISE()
    
    IF (ERROR) THEN
        !WRITE(*,*)'ERROR CALLED - NO HARDENING'
        CALL UPDATE_PARAMETERS()
        RETURN
    ELSE
        CALL HARDENING (FIJ, SIJ, SIJ_NEXT)
        GEFFECTIVE = SECINVDEV(STRESS_TENSOR_NEXT-STRESS_TENSOR_CURRENT)/SECINVDEV(DEPSILON)
        CALL UPDATE_PARAMETERS() 
    END IF

    CONTAINS

    SUBROUTINE OPTIMISE ()
        
        !DSIJ                : INCREMENTAL STRESS RATIO
        !XIJE                : CONJUGATE POINT BASED ON ELASTIC RESPONSE
        !XIJP                : CONJUGATE POINT BASED ON PLASTIC RESPONSE
        !XIJA                : CONJUGATE POINT BASED ON ASSUMPTION
        !XIJC                : SOLVED STRESS TENSOR NEXT
        !SOLN_STRESS_NEXT    : CONTAINS SOLVED STRESS TENSOR NEXT
        !SOLN_SIJ_NEXT        : CONTAINS SOLVED SIJ_NEXT
        
        TYPE(TENSOR_2R)::DSIJ, XIJE, XIJP, XIJA, XIJC, SOLN_STRESS_NEXT, SOLN_SIJ_NEXT
        INTEGER::I
        REAL(KIND=DBL)::THETA_MIN, THETA_CURRENT,SOLN_TANGENT(3,3)
        LOGICAL::FLAG

        FLAG = .FALSE.
        
        ERROR = .FALSE.
        
        !DSIJ = NORMALISE(DEVIATORIC(DEPSILON))
        DSIJ = PREF/HYDROSTATIC_PRESSURE * TENSORMULTIPLY(CONSTITUTIVE_TENSOR,DEVIATORIC(DEPSILON))
        DSIJ = NORMALISE(DSIJ)
        XIJE = SOLVECONJUGATE(SIJ,DSIJ)
        XIJP = RADIUS * NORMALISE(DEVIATORIC(DEPSILON))
                      
        THETA_MIN = 2D0

        DO I = 1, NDIVISION
          
          XIJA = REAL(I,KIND=DBL)/REAL(NDIVISION,KIND=DBL) * (XIJE-XIJP) + XIJP
          XIJA = RADIUS * NORMALISE(XIJA)
          DSIJ = NORMALISE(XIJA-SIJ)

          IF (NORM(SIJ).LE.TOLERANCE.AND.KEY.EQ.0) THEN        ! LOADING POTENTIAL: DEVIATORIC COMPONENT
            FIJ = NORMALISE(DEPSILON)
          ELSE
            FIJ = NORMALISE(SIJ-CENTRE(SURFACE))
          END IF
          
          IF ((FIJ.DDOT.DSIJ).LT.0D0) THEN
            ! UNLOADING/RELOADING (OR INTERPRETED AS SUCH)
            SIJ_NEXT = CALCULATE_STRESS(GIJ=NORMALISE(XIJA), FIJ=NORMALISE(XIJA), SIJ=SIJ, ERROR = ERROR, OPTION = 1) 
          ELSE      
            SIJ_NEXT = CALCULATE_STRESS(GIJ=NORMALISE(XIJA), FIJ=NORMALISE(XIJA), SIJ=SIJ, ERROR = ERROR)
          END IF
          
          IF(ERROR) THEN
            EXIT
          END IF
        
          DSIJ = SIJ_NEXT-SIJ
                   
          IF(NORM(DSIJ).LE. TOLERANCE) THEN
            !WRITE(*,*)'FIRSTIF'
            SOLN_SIJ_NEXT = SIJ_NEXT
            SOLN_STRESS_NEXT = STRESS_TENSOR_NEXT
            SOLN_TANGENT = ELPL_TANGENT_NEXT
            FLAG = .TRUE.
            EXIT
          END IF
          
          DSIJ = NORMALISE(DSIJ)
          XIJC = SOLVECONJUGATE(SIJ,DSIJ)

         THETA_CURRENT = 1.0D0-(NORMALISE(XIJC).DDOT.NORMALISE(XIJA))
                 
         IF(THETA_CURRENT.LT.THETA_MIN.OR.DABS(THETA_MIN).LE.TOLERANCE) THEN
           !WRITE(*,*)'SECONDIF'
           THETA_MIN = THETA_CURRENT
           SOLN_STRESS_NEXT = STRESS_TENSOR_NEXT
           SOLN_SIJ_NEXT = SIJ_NEXT
           SOLN_TANGENT = ELPL_TANGENT_NEXT
           FLAG = .TRUE.
         END IF

       END DO
     
       IF(FLAG) THEN
         STRESS_TENSOR_NEXT = SOLN_STRESS_NEXT
         SIJ_NEXT = SOLN_SIJ_NEXT
         ELPL_TANGENT_NEXT = SOLN_TANGENT
       ELSE
         STRESS_TENSOR_NEXT=STRESS_TENSOR_CURRENT
         ELPL_TANGENT_NEXT = ELPL_TANGENT_CURRENT
         ERROR=.TRUE.
         RETURN
       END IF
      
    END SUBROUTINE OPTIMISE

    SUBROUTINE UPDATE_PARAMETERS()
        CALL UPDATE_VOID_RATIO()
        STRAIN_CURRENT=STRAIN_NEXT
        CALL SENDDATA(DATAOUT=STRESS_CURRENT,DSIZE=3,CHOICE=1)
        CALL SENDDATA(DATAOUT=MODEL_PARAMETER,DSIZE=15,CHOICE=2)
        CALL SENDDATA(DATAOUT=HARDENING_PARAMETER_REAL,DSIZE=7*NSURFACE+3,CHOICE=7)
        CALL SENDDATA(DATAOUT=HARDENING_PARAMETER_INT,DSIZE=2)
        CALL SENDDATA(DATAOUT=TANGENT,DSIZE=3)
    END SUBROUTINE UPDATE_PARAMETERS
   
    END SUBROUTINE MODEL_2D

!******************** START MODULE_GLOBAL_DATA *******************************************

SUBROUTINE RECEIVEDATA_REAL(DATAIN, DSIZE, CHOICE)
        INTEGER,INTENT(IN)::DSIZE,CHOICE
        REAL(KIND=DBL),DIMENSION(DSIZE),INTENT(IN)::DATAIN
        REAL(KIND=DBL),DIMENSION(MAXSIZE,MAXSIZE)::TEMP !MAXSIZE IS 2, DEFINED IN THE MODULE TENSOR_HANDLE
        INTEGER::I,J,K,L
        SELECT CASE (CHOICE)
          CASE(1)
            TEMP=0D0
            TEMP(1,1)=DATAIN(1) ! XX
            TEMP(2,2)=DATAIN(2) ! YY
            TEMP(1,2)=DATAIN(3) ! XY
            TEMP(2,1)=TEMP(1,2) ! YX=XY
            STRESS_TENSOR_CURRENT=CREATETENSOR(TEMP,SIZE_PROB) ! SIZE_PROB IS DEFINED IN THE MODULE KIND_DBL
          CASE(2)
            TEMP=0D0
            TEMP(1,1)=DATAIN(1) ! XX
            TEMP(2,2)=DATAIN(2) ! YY
            TEMP(1,2)=DATAIN(3) ! XY
            TEMP(2,1)=TEMP(1,2) ! YX=XY
            STRAIN_TENSOR_CURRENT=CREATETENSOR(TEMP,SIZE_PROB) ! SIZE_PROB IS DEFINED IN THE MODULE KIND_DBL
          CASE(3)
            TEMP=0D0
            TEMP(1,1)=DATAIN(1) ! XX
            TEMP(2,2)=DATAIN(2) ! YY
            TEMP(1,2)=DATAIN(3) ! XY
            TEMP(2,1)=TEMP(1,2) ! YX=XY
            STRAIN_TENSOR_NEXT=CREATETENSOR(TEMP,SIZE_PROB) ! SIZE_PROB IS DEFINED IN THE MODULE KIND_DBL         
          CASE(4)
            VOID_RATIO                =DATAIN(1)
            SHEAR_CONST_A             =DATAIN(2)
            SHEAR_CONST_N             =DATAIN(3)
            POISSON                   =DATAIN(4)
            ETAMAX_CONST_A1           =DATAIN(5)
            ETAMAX_CONST_B1           =DATAIN(6)
            GNMAX_CONST_A2            =DATAIN(7)
            GNMAX_CONST_B2            =DATAIN(8)
            GNMIN_CONST_A3            =DATAIN(9)
            GNMIN_CONST_B3            =DATAIN(10)
            DEGRADATION               =DATAIN(11)
            DILATANCY_COEFF_MU_MONO   =DATAIN(12)
            DILATANCY_COEFF_MU_CYCLIC =DATAIN(13)
            DILATANCY_STRAIN_SC       =DATAIN(14)
            CRITICAL_STRESS_RATIO_M   =DATAIN(15)
            PATM                      =DATAIN(16)
          CASE(5)
            STEADY_STATE_DATA=0D0
            STEADY_STATE_DATA(:,2)=DATAIN    ! ASSIGNS VOID_RATIOS
          CASE(6)
            STEADY_STATE_DATA(:,1)=DATAIN    ! ASSIGNS PRESSURE
          CASE (7)       
            HYDROSTATIC_DATA=0D0
            HYDROSTATIC_DATA(:,2)=DATAIN     ! ASSIGNS VOID_RATIOS
          CASE (8)
            HYDROSTATIC_DATA(:,1)=DATAIN     ! ASSIGNS PRESSURE
          CASE(9)
            DO I=1,NSURFACE
                TEMP= 0D0
                TEMP(1,1) = DATAIN((I-1)*3+1)
                TEMP(2,2) = DATAIN((I-1)*3+2)
                TEMP(1,2) = DATAIN((I-1)*3+3)
                TEMP(2,1) = TEMP(1,2)
                CENTRE(I) = CREATETENSOR(TEMP, SIZE_PROB)
            END DO
            DO I=1,NSURFACE
                TEMP = 0D0
                TEMP(1,1) = DATAIN(3*NSURFACE+(I-1)*3+1)
                TEMP(2,2) = DATAIN(3*NSURFACE+(I-1)*3+2)
                TEMP(1,2) = DATAIN(3*NSURFACE+(I-1)*3+3)
                TEMP(2,1) = TEMP(1,2)
                REVERSAL(I) = CREATETENSOR(TEMP, SIZE_PROB)
            END DO          
            DO I=1,NSURFACE
                GAMMA(I) = DATAIN(6*NSURFACE+I)
            END DO
            MAX_STRAIN     = DATAIN(7*NSURFACE+1)    
            GEFFECTIVE     = DATAIN(7*NSURFACE+2)
            INITIAL_STRAIN = DATAIN(7*NSURFACE+3)
          END SELECT
    END SUBROUTINE RECEIVEDATA_REAL

    SUBROUTINE RECEIVEDATA_INT(DATAIN,DSIZE)
        INTEGER,INTENT(IN)::DSIZE
        INTEGER,DIMENSION(DSIZE),INTENT(IN)::DATAIN
        SURFACE = DATAIN(1)
        KEY     = DATAIN(2) 
    END SUBROUTINE RECEIVEDATA_INT
    
    SUBROUTINE RECEIVEDATA_REAL_A(DATAIN,DSIZE)
        INTEGER,INTENT(IN)::DSIZE
        REAL(KIND=DBL),DIMENSION(DSIZE,DSIZE),INTENT(IN)::DATAIN
        ELPL_TANGENT_CURRENT=DATAIN 
    END SUBROUTINE RECEIVEDATA_REAL_A
                  
    SUBROUTINE SENDDATA_REAL (DATAOUT,DSIZE,CHOICE)
        INTEGER,INTENT(IN)::DSIZE,CHOICE
        REAL(KIND=DBL),DIMENSION(DSIZE),INTENT(OUT)::DATAOUT
        REAL(KIND=DBL),DIMENSION(MAXSIZE,MAXSIZE)::TEMP !MAXSIZE IS 3, DEFINED IN THE MODULE TENSOR_HANDLE
        INTEGER::I,J,K,L
        SELECT CASE (CHOICE)
          CASE(1)
            TEMP=0D0
            CALL GETTENSOR(STRESS_TENSOR_NEXT,TEMP,SIZE_PROB) ! SIZE_PROB IS DEFINED IN THE MODULE KIND_DBL
            DATAOUT(1)=TEMP(1,1) ! XX
            DATAOUT(2)=TEMP(2,2) ! YY
            DATAOUT(3)=TEMP(1,2) ! XY
          CASE(2)
            DATAOUT(1)  =VOID_RATIO
            DATAOUT(2)  =SHEAR_CONST_A
            DATAOUT(3)  =SHEAR_CONST_N
            DATAOUT(4)  =POISSON
            DATAOUT(5)  =ETAMAX_CONST_A1
            DATAOUT(6)  =ETAMAX_CONST_B1
            DATAOUT(7)  =GNMAX_CONST_A2
            DATAOUT(8)  =GNMAX_CONST_B2
            DATAOUT(9)  =GNMIN_CONST_A3
            DATAOUT(10) =GNMIN_CONST_B3
            DATAOUT(11) =DEGRADATION
            DATAOUT(12) =DILATANCY_COEFF_MU_MONO
            DATAOUT(13) =DILATANCY_COEFF_MU_CYCLIC
            DATAOUT(14) =DILATANCY_STRAIN_SC
            DATAOUT(15) =CRITICAL_STRESS_RATIO_M
            DATAOUT(16) =PATM
           CASE(3)
            DATAOUT=STEADY_STATE_DATA(:,2)
          CASE(4)
            DATAOUT=STEADY_STATE_DATA(:,1)
          CASE(5)       
            DATAOUT=HYDROSTATIC_DATA(:,2)
          CASE(6)
            DATAOUT=HYDROSTATIC_DATA(:,1)
          CASE(7)
            DO I=1,NSURFACE
              TEMP=0D0
              CALL GETTENSOR(CENTRE(I),TEMP,SIZE_PROB)
              DATAOUT((I-1)*3+1) = TEMP(1,1)
              DATAOUT((I-1)*3+2) = TEMP(2,2)
              DATAOUT((I-1)*3+3) = TEMP(1,2)
            END DO
            DO I=1,NSURFACE
              TEMP=0D0
              CALL GETTENSOR(REVERSAL(I),TEMP,SIZE_PROB)
              DATAOUT(3*NSURFACE+(I-1)*3+1) = TEMP(1,1)
              DATAOUT(3*NSURFACE+(I-1)*3+2) = TEMP(2,2)
              DATAOUT(3*NSURFACE+(I-1)*3+3) = TEMP(1,2)
            END DO
            DO I=1,NSURFACE
              DATAOUT(6*NSURFACE+I)=GAMMA(I)
            END DO
            DATAOUT(7*NSURFACE+1)=MAX_STRAIN
            DATAOUT(7*NSURFACE+2)=GEFFECTIVE
            DATAOUT(7*NSURFACE+3)=INITIAL_STRAIN
         END SELECT
    END SUBROUTINE SENDDATA_REAL

    SUBROUTINE SENDDATA_INT(DATAOUT,DSIZE)
        INTEGER,INTENT(IN)::DSIZE
        INTEGER,DIMENSION(DSIZE),INTENT(OUT)::DATAOUT
        DATAOUT(1)=SURFACE
        DATAOUT(2)=KEY    
    END SUBROUTINE SENDDATA_INT
    
    SUBROUTINE SENDDATA_REAL_A(DATAOUT,DSIZE)
        INTEGER,INTENT(IN)::DSIZE
        REAL(KIND=DBL),DIMENSION(DSIZE,DSIZE),INTENT(OUT)::DATAOUT
        DATAOUT=ELPL_TANGENT_NEXT
    END SUBROUTINE SENDDATA_REAL_A
             
    SUBROUTINE CALCULATE_HYDROSTATIC_PRESSURE ()
        HYDROSTATIC_PRESSURE=FIRSTINVARIANT(STRESS_TENSOR_CURRENT)/REAL(SIZE_PROB,KIND=DBL)
    END SUBROUTINE CALCULATE_HYDROSTATIC_PRESSURE
        
    SUBROUTINE CALCULATE_ELASTIC_SHEAR_MODULUS ()
        REAL(KIND=DBL)::N, FACTOR
        CALL CALCULATE_HYDROSTATIC_PRESSURE()
        N = SHEAR_CONST_N
        IF(MAX_STRAIN.GE.0.02D0)THEN
          FACTOR=MIN(MAX_STRAIN/0.05D0,1.0D0)
          N = N+(0.85D0-N)*FACTOR
        END IF
        ELASTIC_SHEAR_MODULUS=SHEAR_CONST_A*PATM*(2.17-VOID_RATIO)*(2.17-VOID_RATIO)/(1+VOID_RATIO)* &
        & ((HYDROSTATIC_PRESSURE/PATM)**N)
    END SUBROUTINE CALCULATE_ELASTIC_SHEAR_MODULUS
       
    SUBROUTINE CALCULATE_STEADY_STATE_VOID_RATIO ()
        INTEGER::I
        REAL(KIND=DBL)::ELIM
        ELIM = HYDROSTATIC_DATA(1,2) - 0.001D0
        CALL CALCULATE_HYDROSTATIC_PRESSURE()
        IF(HYDROSTATIC_PRESSURE.LT.STEADY_STATE_DATA(1,1)) THEN
            !STEADY_STATE_VOID_RATIO = STEADY_STATE_DATA(1,2)
            STEADY_STATE_VOID_RATIO = STEADY_STATE_DATA(1,2)+(STEADY_STATE_DATA(1,1)-HYDROSTATIC_PRESSURE)* &
                (ELIM-STEADY_STATE_DATA(1,2))/STEADY_STATE_DATA(1,1)
        ELSEIF (HYDROSTATIC_PRESSURE.GT.STEADY_STATE_DATA(10,1)) THEN
            STEADY_STATE_VOID_RATIO = STEADY_STATE_DATA(10,2)
        ELSE
          DO I=1,9
            IF(HYDROSTATIC_PRESSURE.GE.STEADY_STATE_DATA(I,1).AND.HYDROSTATIC_PRESSURE.LE.STEADY_STATE_DATA(I+1,1)) THEN
              STEADY_STATE_VOID_RATIO = STEADY_STATE_DATA(I,2)+ &
              (HYDROSTATIC_PRESSURE - STEADY_STATE_DATA(I,1))/(STEADY_STATE_DATA(I+1,1)-STEADY_STATE_DATA(I,1))* &
              (STEADY_STATE_DATA(I+1,2)-STEADY_STATE_DATA(I,2))
              EXIT
            END IF
          END DO
        END IF
    END SUBROUTINE CALCULATE_STEADY_STATE_VOID_RATIO

    SUBROUTINE CALCULATE_HYDROSTATIC_VOID_RATIO ()
        INTEGER::I      
        CALL CALCULATE_HYDROSTATIC_PRESSURE()
        IF(HYDROSTATIC_PRESSURE.LT.HYDROSTATIC_DATA(1,1)) THEN
            HYDROSTATIC_VOID_RATIO = HYDROSTATIC_DATA(1,2)
        ELSEIF (HYDROSTATIC_PRESSURE.GT.HYDROSTATIC_DATA(10,1)) THEN
            HYDROSTATIC_VOID_RATIO = HYDROSTATIC_DATA(10,2)
        ELSE
          DO I=1,9
            IF(HYDROSTATIC_PRESSURE.GE.HYDROSTATIC_DATA(I,1).AND.HYDROSTATIC_PRESSURE.LE.HYDROSTATIC_DATA(I+1,1)) THEN
              HYDROSTATIC_VOID_RATIO = HYDROSTATIC_DATA(I,2)+ &
              (HYDROSTATIC_PRESSURE - HYDROSTATIC_DATA(I,1))/(HYDROSTATIC_DATA(I+1,1)-HYDROSTATIC_DATA(I,1))* &
              (HYDROSTATIC_DATA(I+1,2)-HYDROSTATIC_DATA(I,2))
              EXIT
            END IF
          END DO
        END IF
    END SUBROUTINE CALCULATE_HYDROSTATIC_VOID_RATIO
       
    SUBROUTINE CALCULATE_STATE_INDEX ()
        CALL CALCULATE_STEADY_STATE_VOID_RATIO()
        CALL CALCULATE_HYDROSTATIC_VOID_RATIO()
        STATE_INDEX=(HYDROSTATIC_VOID_RATIO-VOID_RATIO)/(HYDROSTATIC_VOID_RATIO-STEADY_STATE_VOID_RATIO)
    END SUBROUTINE CALCULATE_STATE_INDEX
    
    SUBROUTINE CALCULATE_ETAMAX ()
        CALL CALCULATE_STATE_INDEX()
        ETAMAX=ETAMAX_CONST_A1+ETAMAX_CONST_B1*STATE_INDEX
    END SUBROUTINE CALCULATE_ETAMAX
   
    SUBROUTINE CALCULATE_GNMAX ()
        CALL CALCULATE_STATE_INDEX()
        GNMAX=GNMAX_CONST_A2+GNMAX_CONST_B2*STATE_INDEX
    END SUBROUTINE CALCULATE_GNMAX

    SUBROUTINE CALCULATE_GNMIN ()
        CALL CALCULATE_STATE_INDEX()
        GNMIN=GNMIN_CONST_A3+GNMIN_CONST_B3*STATE_INDEX
    END SUBROUTINE CALCULATE_GNMIN
    
    SUBROUTINE CALCULATE_DELTA_VOLUMETRIC_STRAIN ()
        DELTA_VOLUMETRIC_STRAIN=FIRSTINVARIANT(STRAIN_TENSOR_NEXT)-FIRSTINVARIANT(STRAIN_TENSOR_CURRENT)
    END SUBROUTINE CALCULATE_DELTA_VOLUMETRIC_STRAIN
  
    SUBROUTINE UPDATE_VOID_RATIO ()
        CALL CALCULATE_DELTA_VOLUMETRIC_STRAIN
        VOID_RATIO=VOID_RATIO*(1D0-DELTA_VOLUMETRIC_STRAIN)
    END SUBROUTINE UPDATE_VOID_RATIO

    SUBROUTINE CALCULATE_RADIUS ()
        CALL CALCULATE_ETAMAX
        RADIUS=2D0**0.5D0*ETAMAX*PREF
    END SUBROUTINE CALCULATE_RADIUS  

    SUBROUTINE GET_STATE_DATA(DATASTORE, SIZE)
        INTEGER,INTENT(IN)::SIZE
        REAL(KIND=DBL),DIMENSION(SIZE),INTENT(OUT)::DATASTORE
        DATASTORE(1)=STEADY_STATE_VOID_RATIO
        DATASTORE(2)=HYDROSTATIC_VOID_RATIO
        DATASTORE(3)=STATE_INDEX
    END SUBROUTINE GET_STATE_DATA 

    TYPE(TENSOR_2R) FUNCTION SOLVECONJUGATE_TYPE1 (SIJ, NIJ)      
        TYPE(TENSOR_2R),INTENT(IN)::SIJ,NIJ
        REAL(KIND=DBL)::A, B, C, DETERMINANT, ROOT
        IF (NORM(NIJ).LE.TOLERANCE) THEN
            SOLVECONJUGATE_TYPE1 = SIJ
            RETURN
        END IF
        CALL CALCULATE_RADIUS
        ! PARAMETERS OF THE QUADRATIC
        A = 1.D0
        B = 2.D0*(SIJ.DDOT.NIJ)
        C = (SIJ.DDOT.SIJ)-RADIUS**2
        DETERMINANT=B*B-4.D0*A*C
        IF (DETERMINANT.LT.0D0) THEN
            WRITE(*,*)'STRESS DILATANCY MODEL::SOLVECONJUGATE, DETERMINANT < 0'
            ROOT=0D0
            STOP
        ELSE
            ROOT=(-B+DETERMINANT**0.5D0)/(2.D0*A)
        END IF    
        SOLVECONJUGATE_TYPE1=ROOT*NIJ+SIJ
    END FUNCTION SOLVECONJUGATE_TYPE1
    
    SUBROUTINE SETINITIALTANGENT()
    
        REAL(KIND=DBL)::MU,G
        
        EL_TANGENT = 0D0
        MU = POISSON
        G  = ELASTIC_SHEAR_MODULUS
        
        !-----------------------------------------------------------
        ! THE MATRIX IS IN TERMS OF ENGINEERING STRAIN
        !-----------------------------------------------------------
        
        EL_TANGENT(1,1)=2D0*G/(1D0-2D0*MU)*(1D0-MU)
        EL_TANGENT(1,2)=2D0*G/(1D0-2D0*MU)*MU
                
        EL_TANGENT(2,1)=2D0*G/(1D0-2D0*MU)*MU
        EL_TANGENT(2,2)=2D0*G/(1D0-2D0*MU)*(1D0-MU)
                    
        EL_TANGENT(3,3)=G
        
    END SUBROUTINE SETINITIALTANGENT

    SUBROUTINE GET_REF_PRESSURES()

        !PREF = PATM/100.0D0  
        PREF = 1.0D0
        PMIN = PATM/1000000.0D0

    END SUBROUTINE GET_REF_PRESSURES
    
    REAL(KIND=DBL) FUNCTION CALCULATE_EQV_STRAIN(ETA)
    
        REAL(KIND=DBL),INTENT(IN)::ETA
        REAL(KIND=DBL)::STRAIN_INITIAL,STRAIN_FINAL,GN
        INTEGER::COUNTER
        
        IF(ETA.LE.TOLERANCE) THEN
            CALCULATE_EQV_STRAIN=0D0
            RETURN
        END IF
        
        STRAIN_INITIAL = 0D0
        STRAIN_FINAL   = 0D0
        COUNTER        = 0
        
        DO
        
            GN = (GNMAX-GNMIN)*DEXP(-DEGRADATION*STRAIN_INITIAL/0.01D0)+GNMIN
            STRAIN_FINAL = ETA/(GN*(1D0-ETA/ETAMAX))
            COUNTER = COUNTER + 1
            
            IF(DABS(STRAIN_FINAL-STRAIN_INITIAL).LE.1D-6) EXIT
            
            IF(COUNTER.GT.100) EXIT
            
            STRAIN_INITIAL = STRAIN_FINAL
            
       END DO
                                                             
       CALCULATE_EQV_STRAIN = 0.5D0*(STRAIN_INITIAL+STRAIN_FINAL)
    
    END FUNCTION CALCULATE_EQV_STRAIN

!******************** END MODULE_GLOBAL_DATA *******************************************

!******************** START MODULE_TENSOR_HANDLE ***************************************

       TYPE(TENSOR_2R) FUNCTION TSUM (DATA1, DATA2)
        TYPE(TENSOR_2R),INTENT(IN)::DATA1,DATA2
        IF (DATA1%TSIZE.NE.DATA2%TSIZE)THEN
              WRITE(*,*),'STRESS DILATANCY MODEL::TENSOR_HANDLE::SUM, TENSORS ARE OF DIFFERENT SIZES'
              STOP
          ELSE
              TSUM%TSIZE=DATA1%TSIZE
            TSUM%COMPONENT=DATA1%COMPONENT+DATA2%COMPONENT
          END IF
    END FUNCTION TSUM
    
    TYPE (TENSOR_2R) FUNCTION SUBTRACT (DATA1, DATA2)
        TYPE(TENSOR_2R),INTENT(IN)::DATA1,DATA2
        IF (DATA1%TSIZE.NE.DATA2%TSIZE) THEN
               WRITE(*,*),'STRESS DILATANCY MODEL::TENSOR_HANDLE::SUBTRACT, TENSORS ARE OF DIFFERENT SIZES'
               STOP
        ELSE
               SUBTRACT%TSIZE=DATA1%TSIZE
            SUBTRACT%COMPONENT=DATA1%COMPONENT-DATA2%COMPONENT
        END IF
    END FUNCTION SUBTRACT

    REAL (KIND=DBL) FUNCTION CONTRACT(DATA1,DATA2)
        TYPE(TENSOR_2R),INTENT(IN)::DATA1,DATA2
        IF (DATA1%TSIZE.NE.DATA2%TSIZE)THEN
              WRITE(*,*),'STRESS DILATANCY MODEL::TENSOR_HANDLE::CONTRACT, TENSORS ARE OF DIFFERENT SIZE'     
            STOP
        ELSE
            CONTRACT = SUM(DATA1%COMPONENT*DATA2%COMPONENT)
        END IF
    END FUNCTION CONTRACT  
    
    TYPE (TENSOR_2R) FUNCTION MULTIPLY_DOUBLE_REAL (DATA1, DATA2)
        REAL(KIND=DBL),INTENT(IN)::DATA1
        TYPE(TENSOR_2R),INTENT(IN)::DATA2
        MULTIPLY_DOUBLE_REAL%TSIZE=DATA2%TSIZE
        MULTIPLY_DOUBLE_REAL%COMPONENT=DATA1*DATA2%COMPONENT
    END FUNCTION MULTIPLY_DOUBLE_REAL

    TYPE (TENSOR_2R) FUNCTION MULTIPLY_SINGLE_REAL (DATA1, DATA2)
        REAL,INTENT(IN)::DATA1
        TYPE(TENSOR_2R),INTENT(IN)::DATA2
        MULTIPLY_SINGLE_REAL%TSIZE=DATA2%TSIZE
        MULTIPLY_SINGLE_REAL%COMPONENT=DATA1*DATA2%COMPONENT
    END FUNCTION MULTIPLY_SINGLE_REAL
   
    TYPE (TENSOR_2R) FUNCTION CREATETENSOR (TENSORDATA,DSIZE)
        INTEGER,INTENT(IN)::DSIZE
        REAL(KIND=DBL),DIMENSION(DSIZE,DSIZE),INTENT(IN)::TENSORDATA
        IF (DSIZE.GT.MAXSIZE) THEN
              WRITE(*,*),'STRESS DILATANCY MODEL::TENSOR_HANDLE::CREATETENSOR, EXCEEDS MAXIMUM DIMENSIONS OF 2'
              STOP
        ELSE
              CREATETENSOR%TSIZE=DSIZE
              CREATETENSOR%COMPONENT=TENSORDATA
        END IF
    END FUNCTION CREATETENSOR

    SUBROUTINE GETTENSOR (TENSORDATA,ARRAYDATA,DSIZE)
        TYPE (TENSOR_2R),INTENT(IN)::TENSORDATA
        INTEGER,INTENT(IN)::DSIZE
        REAL(KIND=DBL),DIMENSION(DSIZE,DSIZE),INTENT(OUT)::ARRAYDATA
        ARRAYDATA=TENSORDATA%COMPONENT
   END SUBROUTINE GETTENSOR

    TYPE (TENSOR_2R) FUNCTION NORMALISE(TENSORDATA)
         TYPE(TENSOR_2R),INTENT(IN)::TENSORDATA
        REAL(KIND=DBL)::NORM
        NORM=CONTRACT(TENSORDATA,TENSORDATA)**0.5D0
        IF (NORM.LE.TOLERANCE) THEN
            ! DON'T DIVIDE BY ZERO!!!!
            NORMALISE%TSIZE = TENSORDATA%TSIZE
            NORMALISE%COMPONENT = 0.0D0
        ELSE
            NORMALISE=1D0/NORM*TENSORDATA
        END IF
    END FUNCTION NORMALISE

    REAL(KIND=DBL) FUNCTION FIRSTINVARIANT (TENSORDATA)
         TYPE(TENSOR_2R),INTENT(IN)::TENSORDATA
        INTEGER::I
        FIRSTINVARIANT=0D0
        DO I=1,TENSORDATA%TSIZE
          FIRSTINVARIANT=TENSORDATA%COMPONENT(I,I)+FIRSTINVARIANT
        END DO
    END FUNCTION FIRSTINVARIANT

    TYPE (TENSOR_2R) FUNCTION HYDROSTATIC (TENSORDATA)
         TYPE(TENSOR_2R),INTENT(IN):: TENSORDATA
        INTEGER::I
        REAL(KIND=DBL)::INV1
        INV1=FIRSTINVARIANT(TENSORDATA)/REAL(TENSORDATA%TSIZE)
        HYDROSTATIC%TSIZE=TENSORDATA%TSIZE
        HYDROSTATIC%COMPONENT=0D0
        DO I=1,HYDROSTATIC%TSIZE
              HYDROSTATIC%COMPONENT(I,I)= INV1
        END DO
    END FUNCTION HYDROSTATIC

    TYPE (TENSOR_2R) FUNCTION DEVIATORIC (TENSORDATA)
         TYPE(TENSOR_2R),INTENT(IN):: TENSORDATA
        INTEGER::I 
        REAL(KIND=DBL)::INV1   
        INV1=FIRSTINVARIANT(TENSORDATA)/REAL(TENSORDATA%TSIZE,KIND=DBL)
        DEVIATORIC%TSIZE=TENSORDATA%TSIZE
        DEVIATORIC%COMPONENT=TENSORDATA%COMPONENT
        DO I=1,DEVIATORIC%TSIZE
          DEVIATORIC%COMPONENT(I,I)=DEVIATORIC%COMPONENT(I,I)-INV1
        END DO
    END FUNCTION DEVIATORIC

    REAL(KIND=DBL) FUNCTION SECINVDEV (TENSORDATA)
        TYPE (TENSOR_2R), INTENT(IN):: TENSORDATA
        SECINVDEV=((DEVIATORIC(TENSORDATA).DDOT.DEVIATORIC (TENSORDATA))/2D0)**0.5D0
    END FUNCTION SECINVDEV

    TYPE (TENSOR_2R) FUNCTION TENSORMULTIPLY222 (DATA1,DATA2,DUMMY1,DUMMY2)
        TYPE (TENSOR_2R),INTENT(IN)::DATA1,DATA2
        ! DUMMY1 IS DUMMY INDEX OF THE TENSOR DATA1
        ! DUMMY2 IS DUMMY INDEX OF THE TENSOR DATA2
        INTEGER,INTENT(IN)::DUMMY1,DUMMY2
        INTEGER::I,J,K
        IF (DATA1%TSIZE.NE.DATA2%TSIZE) THEN
              WRITE(*,*),'STRESS DILATANCY MODEL::TENSOR_HANDLE::TENSORMULTIPLY222, TENSORS ARE OF DIFFERENT SIZES'
              STOP
        END IF
        TENSORMULTIPLY222%TSIZE=DATA1%TSIZE
        DO I=1,DATA1%TSIZE
        DO J=1,DATA1%TSIZE
        !FORALL(I=1:DATA1%TSIZE,J=1:DATA1%TSIZE) TENSORMULTIPLY222%COMPONENT(I,J)=0D0
            TENSORMULTIPLY222%COMPONENT(I,J)=0D0
        END DO 
        END DO
        SELECT CASE (DUMMY1)
          CASE(1)
              IF (DUMMY2.EQ.1) THEN
                   DO I=1,DATA1%TSIZE
                      DO J=1,DATA1%TSIZE
                          DO K=1,DATA1%TSIZE
                            TENSORMULTIPLY222%COMPONENT(I,J)=DATA1%COMPONENT(K,I)*DATA2%COMPONENT(K,J) + &
                            & TENSORMULTIPLY222%COMPONENT(I,J)
                        END DO
                    END DO
                 END DO
            ELSE
                  DO I=1,DATA1%TSIZE
                      DO J=1,DATA1%TSIZE
                          DO K=1,DATA1%TSIZE
                            TENSORMULTIPLY222%COMPONENT(I,J)=DATA1%COMPONENT(K,I)*DATA2%COMPONENT(J,K) + &
                            & TENSORMULTIPLY222%COMPONENT(I,J)
                        END DO
                    END DO
                END DO
            END IF
          CASE(2)
            IF (DUMMY2.EQ.1) THEN
                DO I=1,DATA1%TSIZE
                      DO J=1,DATA1%TSIZE
                          DO K=1,DATA1%TSIZE
                            TENSORMULTIPLY222%COMPONENT(I,J)=DATA1%COMPONENT(I,K)*DATA2%COMPONENT(K,J) + &
                            & TENSORMULTIPLY222%COMPONENT(I,J)
                        END DO
                     END DO
                END DO
            ELSE
                DO I=1,DATA1%TSIZE
                      DO J=1,DATA1%TSIZE
                          DO K=1,DATA1%TSIZE
                            TENSORMULTIPLY222%COMPONENT(I,J)=DATA1%COMPONENT(I,K)*DATA2%COMPONENT(J,K) + &
                            & TENSORMULTIPLY222%COMPONENT(I,J)
                        END DO
                    END DO
               END DO
            END IF             
        END SELECT
    END FUNCTION TENSORMULTIPLY222

    TYPE(TENSOR_4R) FUNCTION CONSTITUTIVETENSOR (SHEAR_MODULUS, MU, DSIZE)
        REAL(KIND=DBL),INTENT(IN)::SHEAR_MODULUS, MU
        INTEGER,INTENT(IN):: DSIZE
        REAL(KIND=DBL),DIMENSION(DSIZE,DSIZE)::KDELTA !KRONECKER'S DELTA
        INTEGER::I,J,K,L
        KDELTA=0D0
        DO I=1,DSIZE
            DO J=1,DSIZE
                IF (I.EQ.J) THEN
                    KDELTA(I,J)=1D0
                ELSE
                    KDELTA(I,J)=0D0
                END IF
            END DO
        END DO
        CONSTITUTIVETENSOR%TSIZE=DSIZE
        CONSTITUTIVETENSOR%COMPONENT=0D0
        DO I=1,DSIZE
            DO J=1,DSIZE
                DO K=1,DSIZE
                    DO L=1,DSIZE
                        CONSTITUTIVETENSOR%COMPONENT(I,J,K,L)=2D0*SHEAR_MODULUS*MU/(1-2D0*MU)*KDELTA(I,J)*KDELTA(K,L)+&
                        & SHEAR_MODULUS*(KDELTA(I,K)*KDELTA(J,L)+KDELTA(I,L)*KDELTA(J,K))                    
                    END DO
                END DO
             END DO
        END DO
    END FUNCTION CONSTITUTIVETENSOR

    TYPE(TENSOR_2R) FUNCTION TENSORMULTIPLY422(DATA1,DATA2)
        TYPE(TENSOR_4R),INTENT(IN)::DATA1
        TYPE(TENSOR_2R),INTENT(IN)::DATA2
        INTEGER::I,J,K,L
        IF(DATA1%TSIZE.NE.DATA2%TSIZE) THEN
          PRINT *,'STRESS DILATANCY MODEL:: TENSOR_HANDLE::TENSORMULTIPLY422, TENSORS ARE OF DIFFERENT SIZE'
          STOP
        END IF
        TENSORMULTIPLY422%TSIZE=DATA1%TSIZE
        TENSORMULTIPLY422%COMPONENT=0D0
        DO I=1,DATA1%TSIZE
            DO J=1,DATA1%TSIZE
                DO K=1,DATA1%TSIZE
                    DO L=1,DATA1%TSIZE
                        TENSORMULTIPLY422%COMPONENT(I,J)=TENSORMULTIPLY422%COMPONENT(I,J)+ &
                        & DATA1%COMPONENT(I,J,K,L)*DATA2%COMPONENT(K,L)
                     END DO
                END DO
            END DO
        END DO
        
        
     END FUNCTION TENSORMULTIPLY422

     REAL(KIND=DBL) FUNCTION NORM (TENSORDATA)
         TYPE (TENSOR_2R),INTENT(IN)::TENSORDATA
        NORM = (TENSORDATA.DDOT.TENSORDATA)**0.5D0
     END FUNCTION NORM

!******************** END MODULE_TENSOR_HANDLE *****************************************

!******************** START MODULE_PLASTIC_FLOW ****************************************

TYPE(TENSOR_2R) FUNCTION CALCULATE_STRESS(GIJ, FIJ, SIJ, ERROR, OPTION)

        TYPE(TENSOR_2R),INTENT(IN)::GIJ,FIJ,SIJ
        INTEGER, INTENT(IN), OPTIONAL::OPTION
        LOGICAL,INTENT(INOUT)::ERROR

        REAL(KIND=DBL)::C,Q,P,TEMP(MAXSIZE,MAXSIZE),DIR

        TYPE(TENSOR_2R)::LOADING_POTENTIAL,LOADING_POTENTIAL_DEVIATORIC, LOADING_POTENTIAL_HYDROSTATIC, &
        & PLASTIC_POTENTIAL,PLASTIC_POTENTIAL_DEVIATORIC,PLASTIC_POTENTIAL_HYDROSTATIC,DEPSILON

        INTEGER::I,J

        !************************************************************************************************
        ! CALCULATE COAXIALITY AND THE Q/P TERM FOR DILATANCY
        !************************************************************************************************

        IF (NORM(SIJ).LE. TOLERANCE) THEN    ! ASSUME SIJ AND GIJ ARE PARALLEL
            C = 1D0
        ELSE
            C = NORMALISE(SIJ).DDOT.GIJ        ! NOTE: GIJ IS ALREADY NORMALISED
        END IF

        Q=SECINVDEV(STRESS_TENSOR_CURRENT)    ! NOTE: Q = (J2)^0.5
        P=HYDROSTATIC_PRESSURE

        IF(P.LE.PMIN) THEN
          Q = ETAMAX*P
          C = 1D0
        END IF
                  
        !************************************************************************************************
        ! CALCULATE PLASTIC AND LOADING POTENTIALS
        !************************************************************************************************
                           
        PLASTIC_POTENTIAL_DEVIATORIC = REAL((1D0/2D0)**0.5D0,KIND=DBL)* GIJ   

        FORALL(I=1:MAXSIZE,J=1:MAXSIZE) TEMP(I,J) = 0D0
        !TEMP = 0D0

        IF(PRESENT(OPTION)) THEN   
            ! THIS MODIFICATION IS TO ACCOUNT FOR INITIAL ANISOTROPIC STRAIN
            !WRITE(*,*)'OPTION PRESENT'     
            DO I=1,SIZE_PROB
                DO J=1,SIZE_PROB
                    IF(I.EQ.J) THEN
                        TEMP(I,J) = (DILATANCY_AT_REVERSAL-C*Q/P)/REAL(SIZE_PROB,KIND=DBL)
                    END IF
                END DO
            END DO
            PLASTIC_SHEAR_MODULUS = PLASTIC_MODULUS_REVERSAL
        ELSE
            DO I=1,SIZE_PROB
                DO J=1,SIZE_PROB
                    IF(I.EQ.J) TEMP(I,J) = (DILATANCY-C*Q/P)/REAL(SIZE_PROB,KIND=DBL)  
                END DO
            END DO           
            
        END IF         
        PLASTIC_POTENTIAL_HYDROSTATIC = CREATETENSOR(TEMP,SIZE_PROB) 
        PLASTIC_POTENTIAL = PLASTIC_POTENTIAL_DEVIATORIC+PLASTIC_POTENTIAL_HYDROSTATIC

        LOADING_POTENTIAL_DEVIATORIC = REAL((1D0/2D0)**0.5D0,KIND=DBL)* FIJ
        DO I=1,SIZE_PROB
            DO J=1,SIZE_PROB
                IF(I.EQ.J) TEMP(I,J) = -C*Q/P/REAL(SIZE_PROB,KIND=DBL) 
            END DO
        END DO  
        LOADING_POTENTIAL_HYDROSTATIC=CREATETENSOR(TEMP,SIZE_PROB)
        LOADING_POTENTIAL=LOADING_POTENTIAL_DEVIATORIC+LOADING_POTENTIAL_HYDROSTATIC

        !************************************************************************************************
        ! CALCULATE LAMDA
        !************************************************************************************************
        
        DEPSILON=STRAIN_TENSOR_NEXT-STRAIN_TENSOR_CURRENT
       
        LAMDA=(LOADING_POTENTIAL.DDOT.TENSORMULTIPLY(CONSTITUTIVE_TENSOR,DEPSILON))/ &
        & (PLASTIC_SHEAR_MODULUS+(LOADING_POTENTIAL.DDOT.TENSORMULTIPLY(CONSTITUTIVE_TENSOR,PLASTIC_POTENTIAL)))

        IF (LAMDA.LT.0D0) LAMDA = 0.0D0

        PLASTIC_STRAIN_INCREMENT=LAMDA*PLASTIC_POTENTIAL

        !************************************************************************************************
        ! CALCULATE NEW STRESS STATE
        !************************************************************************************************
        
        STRESS_TENSOR_NEXT=STRESS_TENSOR_CURRENT+ &
        & TENSORMULTIPLY(CONSTITUTIVE_TENSOR,(DEPSILON-PLASTIC_STRAIN_INCREMENT))

        !************************************************************************************************
        ! PROJECT DEVIATORIC COMPONENTS ON TO THE STANDARD PLANE
        !************************************************************************************************
       
        P=FIRSTINVARIANT(STRESS_TENSOR_NEXT)/REAL(SIZE_PROB,KIND=DBL)
        
        IF(P.LT.TOLERANCE) THEN
            STRESS_TENSOR_NEXT = STRESS_TENSOR_CURRENT
            ERROR = .TRUE.
            RETURN
        END IF
        
        Q = SECINVDEV(STRESS_TENSOR_NEXT)

        IF(Q.GT.P*ETAMAX) THEN 
            CALL STRESSCORRECTION(DEPSILON)
        ELSE
            CALL SETELASTOPLASTICTANGENT(PLASTIC_POTENTIAL,LOADING_POTENTIAL)
        END IF

        P=FIRSTINVARIANT(STRESS_TENSOR_NEXT)/REAL(SIZE_PROB,KIND=DBL)  

        CALCULATE_STRESS=PREF/P*DEVIATORIC(STRESS_TENSOR_NEXT)  
     
    END FUNCTION CALCULATE_STRESS

    SUBROUTINE STRESSCORRECTION(DEPS)
        TYPE(TENSOR_2R),INTENT(IN)::DEPS
        TYPE(TENSOR_2R)::NIJ
        REAL(KIND=DBL)::TEMP(2,2),Q,DELQ,P,DELP,ETA
        DELQ = GEFFECTIVE*SECINVDEV(DEPS)
        ETA  = SECINVDEV(STRESS_TENSOR_CURRENT)/HYDROSTATIC_PRESSURE
        DELP = DELQ/MAX(ETA,ETAMAX)
        Q    = SECINVDEV(STRESS_TENSOR_CURRENT)+DELQ
        P    = HYDROSTATIC_PRESSURE+DELP
        NIJ  = NORMALISE(DEVIATORIC(STRESS_TENSOR_NEXT))
        TEMP=0D0
        TEMP(1,1)=P
        TEMP(2,2)=P
        STRESS_TENSOR_NEXT = (2D0**0.5D0)*Q*NIJ + CREATETENSOR(TEMP,SIZE_PROB)
    END SUBROUTINE STRESSCORRECTION
    
    SUBROUTINE SETELASTOPLASTICTANGENT(PLASTIC_POTENTIAL,LOADING_POTENTIAL)
        TYPE(TENSOR_2R),INTENT(IN)::PLASTIC_POTENTIAL,LOADING_POTENTIAL
        REAL(KIND=DBL)::TEMP(SIZE_PROB,SIZE_PROB),DELG(3,1),DELF(3,1),DENOMINATOR,CORRECTOR(3,3)
        DELG = 0D0 ! VECTOR FORM OF PLASTIC_POTENTIAL
        DELF = 0D0 ! VECTOR FORM OF LOADING_POTENTIAL
        CALL GETTENSOR(PLASTIC_POTENTIAL,TEMP,SIZE_PROB)
        DELG(1,1)=TEMP(1,1)
        DELG(2,1)=TEMP(2,2)
        DELG(3,1)=2D0*TEMP(1,2)
        CALL GETTENSOR(LOADING_POTENTIAL,TEMP,SIZE_PROB)
        DELF(1,1)=TEMP(1,1)
        DELF(2,1)=TEMP(2,2)
        DELF(3,1)=2D0*TEMP(1,2)
        DENOMINATOR = PLASTIC_SHEAR_MODULUS+(LOADING_POTENTIAL.DDOT.TENSORMULTIPLY(CONSTITUTIVE_TENSOR,PLASTIC_POTENTIAL))
        CORRECTOR = MATMUL(EL_TANGENT,MATMUL(MATMUL(DELG,TRANSPOSE(DELF)),EL_TANGENT))
        CORRECTOR = CORRECTOR/DENOMINATOR
        ELPL_TANGENT_NEXT=EL_TANGENT-CORRECTOR
    END SUBROUTINE SETELASTOPLASTICTANGENT

!******************** END MODULE_PLASTIC_FLOW ******************************************

!******************** START MODULE_HARDENING_RULE **************************************

SUBROUTINE HARDENING(FIJ,SIJ,SIJ_TRIAL)     
        TYPE(TENSOR_2R),INTENT(IN)::FIJ, SIJ, SIJ_TRIAL
        TYPE(TENSOR_2R):: DSIJ
        REAL(KIND=DBL):: DGAMMA
                
        DSIJ    = SIJ_TRIAL-SIJ
                           
        !*******************************************************************************************
        ! SWITCH ON KEY
        !*******************************************************************************************
        
        IF (KEY.EQ.0) THEN
            KEY = 1
        END IF

        !*******************************************************************************************
        ! PLASTIC STRAIN INCREMENT
        !*******************************************************************************************
        
        DGAMMA = 2.0D0*SECINVDEV(STRAIN_TENSOR_NEXT-STRAIN_TENSOR_CURRENT)
        !WRITE(*,*)'DGAMMA',DGAMMA

        !*******************************************************************************************
        ! CASE OF PLASTIC FLOW: SIJ = SIJ_NEXT
        !*******************************************************************************************

        IF(NORM(DSIJ).LT.TOLERANCE) THEN
            GAMMA(SURFACE)=GAMMA(SURFACE)+DGAMMA
            RETURN
        END IF
                
        !*******************************************************************************************
        ! CASE OF LOADING   : DELETE SMALLER SURFACES AND ACTIVATE THE NEXT LARGER SURFACE
        ! CASE OF UNLOADING : CREATE NEW SURFACE
        !*******************************************************************************************
    
        IF((DSIJ.DDOT.FIJ).GE.0D0) THEN
                                                                ! LOADING    
            CALL ACTIVATE_SURFACE(SURFACE,SIJ_TRIAL,DGAMMA)        
            !WRITE(*,*)'HARDENING-LOADING'
                                                                ! POPULATE / UPDATE CURRENT SURFACE
        ELSE
                                                                  ! UNLOADING 
            MAX_STRAIN=MAX(MAX_STRAIN,GAMMA(SURFACE))  
            !WRITE(*,*)'HARDENING-UNLOADING'
            
            SURFACE = SURFACE+1                                    
                                                                ! CREATES ONE ADDITIONAL SURFACE UPTO NSURFACE+1
            IF (SURFACE.LE.NSURFACE) THEN
              REVERSAL(SURFACE) = SIJ
            ELSE
              REVERSAL(SURFACE-1) = SIJ
            END IF
                                                                  ! ASSIGNS THE REVERSAL STRESS
            CALL ACTIVATE_SURFACE(SURFACE,SIJ_TRIAL,DGAMMA)        
                                                                  ! GENERATE NEW SURFACE LIMITED TO NSURFACE           
       END IF

    END SUBROUTINE HARDENING

    RECURSIVE SUBROUTINE ACTIVATE_SURFACE(N,SIGMA,DGAMMA)
        INTEGER,INTENT(IN)::N
        TYPE(TENSOR_2R),INTENT(IN)::SIGMA
        REAL(KIND=DBL),INTENT(IN)::DGAMMA
        TYPE(TENSOR_2R)::XIJ_CP,NIJ
        REAL(KIND=DBL)::DIST,RADIUS
        REAL(KIND=DBL)::TEMP(MAXSIZE,MAXSIZE)
        INTEGER::I,J
        
        IF(N.EQ.1)THEN 
                                                ! TERMINATE THE RECURSION
            SURFACE = N
            REVERSAL(SURFACE) = SIGMA
            GAMMA(SURFACE) = GAMMA(SURFACE)+DGAMMA
            NIJ=NORMALISE(SIGMA)
            XIJ_CP=SOLVECONJUGATE(SIGMA,NIJ)
            CENTRE(SURFACE)=SIGMA-(NORM(SIGMA)/NORM(XIJ_CP))*(XIJ_CP)
            
        ELSE
            
            DIST =  NORM(SIGMA - CENTRE(N-1))
            RADIUS = NORM(REVERSAL(N-1)- CENTRE(N-1))
            
            IF (DIST.GT.RADIUS) THEN            
                                                ! DELETE THE CURRENT SURFACE IF N <= NSURFACE
                                                ! FOR N = NSURFACE+1 THERE IS NO MEMORY
                                                ! MAKE CURRENT SURFACE = N-1
                IF(N.LE.NSURFACE) THEN   
                                               
                    FORALL(I=1:MAXSIZE,J=1:MAXSIZE) TEMP(I,J)=0D0
                    REVERSAL(N)    = CREATETENSOR(TEMP, SIZE_PROB)
                    CENTRE(N)    = CREATETENSOR(TEMP, SIZE_PROB)
                    GAMMA(N)    = INITIAL_STRAIN
                    
                END IF
                
                CALL ACTIVATE_SURFACE(N-1,SIGMA,0D0)
                
            ELSE                                 ! POPULATE THE CURRENT SURFACE

                  IF(N.LE.NSURFACE) THEN
                
                      SURFACE = N
                      CENTRE(SURFACE) = LOCATE_CENTRE(REVERSAL(SURFACE),CENTRE(SURFACE-1),SIGMA)
                      GAMMA(SURFACE)= GAMMA(SURFACE)+DGAMMA

                ELSE
                                                ! DUMMY SURFACE =  NSURFACE + 1 IS RE-ASSIGNED AS
                                                ! SURFACE = NSURFACE; SIJ IS THE STRESS REVERSAL OF
                                                ! THE DUMMY SURFACE
                    SURFACE = N-1
                    CENTRE(SURFACE) = LOCATE_CENTRE(REVERSAL(SURFACE), CENTRE(SURFACE), SIGMA)
                    GAMMA(SURFACE) = GAMMA(SURFACE)+DGAMMA

                END IF
             
            END IF
            
        END IF
        
     END SUBROUTINE ACTIVATE_SURFACE

     TYPE(TENSOR_2R) FUNCTION LOCATE_CENTRE(SIGMA_R, SIGMA_C, SIGMA) RESULT(Z)
         TYPE(TENSOR_2R),INTENT(IN)::SIGMA_R, SIGMA_C, SIGMA
        TYPE(TENSOR_2R)::NIJ
        REAL(KIND=DBL)::BETA,NUMERATOR,DENOMINATOR

        NIJ = NORMALISE(SIGMA_R-SIGMA_C)

        NUMERATOR = (SIGMA_R.DDOT.(SIGMA_R-SIGMA_C))+(SIGMA.DDOT.(SIGMA_C-SIGMA))+(SIGMA_C.DDOT.(SIGMA-SIGMA_R))

        DENOMINATOR = NIJ.DDOT.(SIGMA_R-SIGMA)

         IF(DABS(DENOMINATOR).LE.TINY(1D0)) THEN
            Z = SIGMA_C
            RETURN
         ELSE    
            BETA = 0.5D0*NUMERATOR/((NIJ.DDOT.(SIGMA_R-SIGMA)))
            Z = BETA*(NIJ)+SIGMA_C
         END IF
        
     END FUNCTION LOCATE_CENTRE

!******************** END MODULE_HARDENING_RULE ****************************************

END MODULE ALL_INTERFACES_2D
